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INTRODUCTION 

The  present  study  arose  in  connection  with,  a  math- 
ematical problem  in  the  field  of  diffusion  theory,  and 
was  undertaken  by  the  New  York  University  Group  of  the 
Applied  Mathematics  Panel  at  the  request  of  Division  10  of 
the  N.D.R.C   It  was  carried  out  by  Dr.  B.  Friedman  and 
Dr.  M.  Shiffman  with  the  assistance  of  S.  Isaacson. 

An  aerosol  is  a  suspension  of  liquid  particles  in  air, 
such  as  vapor  or  smoke.   By  mutual  collisions,  particles  will 
coagulate;  by  gravity,  particles  will  settle  on  the  ground 
and  disappear  from  the  aerosol.   The  problem  is  to  determine 
the  number   N(R,T)   of  particles  of  radius  R   at  any  time   T, 
provided  that  the  probability  of  collision  and  that  of  settling 
are  known"  and  that  the  number   N  of  particles  of  radius   R 
at  the  time   T  =0   is  given  as  a  function  of  R. 

The  corresponding  mathematical  problem  is  that  of  solv- 
ing the  initial  value  problem  of  a  non-linear  integro-dif f er- 
ential  equation  of  the  ''  convolution  type;i .   This  equation  sets' 
the  rate  of  change  of  the  number   N(R,T)  of  particles  equal 
to  the  rate  of  creation  of  these  particles  by  coagulation  of 
smaller  particles,  minus  the  rate  at  which  the  particles  of 
radius  R   collide  with  other  particles  and  minus  the  rate  at 

■5C-  These  probability  functions  were  determined  long  ago  by 
Smoluchowski  who  seems  to  have  been  the  first  to  investigate 
the  coagulation  problem.   The  settling  term  was  included  by 
Dr.  Goldman  of  Division  10  of  the  N.D.R.C,  who  derived  the 
integral  equation,  the  solution  of  which  is  the  subject  of  the 
present  memorandum* 
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which  they  settle. 

In  the  present  paper  a  solution  of  the  problem  is 
given  provided  the  initial  distribution  N(R,0)  is  homo- 
geneous, i.e.,  all  the  particles  have  the  same  radius.   The 
procedure  is  not  essentially  different  from  Smoluchowski ' s 
method  and  is  based  on  a  reduction  to  a  system  of  ordinary 
differential  equations. 

In  a  subsequent  memorandum  a  solution  of  the  integro- 
differential  equation  for  non-homogeneous  initial  distributions 
will  be  presented.   Such  an  extension  of  the  theory  is  desirable 
since  in  practice  aerosols  may  not  be  initially  homogeneous. 
The  generalized  method  will  likewise  apply  to  a  larger  class 
of  functional  equations  of  diffusion. 


R.  Courant 
Contractor's  Technical  Representative 
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SOLUTION  OF  A  NON-LINEAR  INTEGRAL  EQUATION  APPEARING  IN 

THE  THEORY  OF  AN  AEROSOL 

1.  The  Integral  Equation* 

The  behavior  of  an  aerosol,  i.  e.,  a  suspension  of 
spherical  liquid  particles  in  air,  is  of  interest  in  many 
practical  applications.   The  particles,  in  the  course  of 
their  random  motion,  collide  with  one  another  and  coagulate 
to  form  larger  particles.   In  addition,  some  particles 
settle  to  the  bottom  of  the  air  container  by  the  action  of 
gravity,  and  so  disappear  from  the  suspension.  A  quantitative 
theory  of  this  action  was  developed  by  Dr.  D*  E«  Goldman, 
who  derived  a  non-linear  integro-diff erential  equation  on  the 
basis  of  kinetic  theory.   This  memorandum  takes  the  integral 
equation  as  a  starting  point  and  offers  a  solution.   A  state- 
ment of  Dr.  Goldman's  result  follows; 

Let  N(R,T)dR   be  the  concentration  of  particles  of 
radii  between  R   and  R+dR   at  time   T.   It  'is  required 
to  find  N(R,T)  given  the  initial  distribution  N(R,0). 
For  any  fixed  R,  the  number  of  particles  will  change  with 
time  because  of  the  following  factors; 

1.  Particles  of  size  R  will  be  destroyed  by  collision 
with  particles  of  any  other  size. 

2.  Particles  of  size  R   will  be  created  by  the  coagula- 
tion of  smaller  particles. 

3.  particles,  of  size  R  will  settle  under  the  influence 
of  gravity. 

-::-  In  a  manuscript  dated  11-29-43,  for  Section  10.2,  NDRC 
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Expressing  these  effects  mathematically,  the  follow- 
ing equation  was  obtained  by  Goldman; 


(1) 


3N(R,T)  _ 


/  f2f(X,R)N(R,T)N(X,T)dX 

2  oc 

K    P  *: 

+  -3   /  0(X,Y)M(X,T)1J(Y,T)  • 

4   u  -   „    ., 

X°+Y,3=  R6 


dR 


-  U(R)II(R,T), 


where      #(X,Y)    and     U(R)      are    the   following   functions; 


0(X,Y)    = 


X    +  DC 


OC 


X' 


Y' 


(X  +  Y) 


U(R)    =    /SR(R    +  OC). 

The  constants   oc,  8     ^-nd  K   are  defined  as  follows; 

'  o 

cC    is    connected   with    the   mean  free    path  of   the    air  molecules 
and    can  be    taken   as 

-4 
CX    =   .083  microns   =   .083  X    lu        cm., 

(3     =  2Pg/9l9h      =    1.27  X    106    P/h, 

K      =  4KT/5>7  =  3  X   10"10    (for   T   =  273°   absolute) 

where    P  is    the   density  of    the    liquid,    g    is    the  gravitational 

acceleration,  77  the    viscosity  of   the    air,    h   the   ratio 

of   volume    to    floor   area   of   the    container,   K  is   Boltzman's 


*  In   the   second   integral   the    range    of  X   is  from   0_to   R 
and  Y  is    determined    by   the   equation  Y3   =  R3    -  X3.      This 
equation  expresses   the    fact    that,    in    coagulation,    the 
volumes    add. 
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constant   and   T   the    absolute    temperature.      (The    values   of 
the    constants    are    given    in   c.g.s.   units.) 

In   e que  t ion    (1),    the    individual  terras    on    the    right 
side   represent,   respectively,    the   effect   of   the  destruction, 
creation   and    settling  of    the    particles. 

It   is    convenient    to   introduce    non-dimensional  vari- 
ables  r,   x,   y  as    follows*      Let  ct  be    the    unit    of    length, 
1/vX'"^     the   unit   of    time,    2CC|8/K      the    unit    of   concentration 
per   unit    radius,    and    2  OC^^S /K     the    unit    of    total   concentra- 
tion.     Nov,   set 

R  =  Otr,  X  =-ocx,  Y  =  <Xy 

^N(R,T)   =  2oc/Sn(r,t)/K:o 

A  simple    integration   of   the    last   equation    justifies    our 
choice    of    the    unit   of   total   concentration, 
hen    (2)    is    substituted    in    (1)   we    obtain 


(3)^  n(r,t) 


/ 
/ 


vo 


f (r,x  )n(r,t )n(x,  t  )dx 


+ 


|     /   f(x,y)nU,t)n(y,t)  §2  dx 


where    now 


-   u(r)n(r,t ), 


(4)        f(x,y)    = 


x+1 
pj 

Ct 

X 


y+1   1 

y   J 


(x+y) 
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(5)  u(r)    =  r(r  +   1)    . 

The  problem  before  us  now  is  to  solve  the  integral 
equation  (3)  given  the  initial  distribution  n(r,0). 

The  equation  (3)  is  a  'reduced'  equation,  which  does 
not  involve  any  physical  constants*.   This  shows  that  if  two 
aerosols  have  the  same  reduced  initial  distribution  n(r,0), 
they  have  similar  behavior* 

2»  Summary  of  the  Results 

In  this  memorandum  the  problem  is  solved  under  the 
assumption  that  the  initial  distribution  in  the  aerosol  is 
homogeneous,  i.e.,  particles  of  only  one  size  are  present. 

Let  the  initial  distribution  have  a  (non-dimensional) 
concentration  n-,   of  particles  of  a  single  radius   r,. 
B'Y  coagulation,  particles  of  twice  the  mass,  or  of  radius 
yj/2r1   arc  formed.   By  further  collisions  and  coagulations, 
particles  of  radii  "^~5v-,f      ^[Tr^,  etc.  .will  appear.  ,  For 
convenience  denote  ^f^^i,      ^J^>Tlt      ^J^rif    etc.,  by 
r0,  r,,  r,,  etc.  respectively. 

Let  n-,(t),  rv,(t),  n~(t),  etc.,  denote  the  concentration 
of  particles,  at  any  time  t,  of  radii   r, ,  r?,  r„etc,  respect- 
ively.  The  particles  of  radius   rg   are  created  rapidly  at 
the  beginning,  while  those  of  radii   r^,  r^,  etc.,  are  created 
very  slowly.   For  practical  purposes  we  can  therefore  set 
n,(t)  =  n4(t)  =....=  0. 

The  solutions  for   n-,(t)   and  ng(t)   are  divided  in- 
to two  cases  depending  on  whether  the  settling  or  collision 
term  is  initially  predominant.   If  the  settling  term  pre- 
dominates over  the  collision  term,  we  obtain  the  formulas; 
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(6)    ni(t) 


n-^l  -  A*) 


(7)    n2(t)  =  ^(t)-!  Ae 


-.  6un  t        [ 


u. 


± 


log 


1  +  A   e  c       -  A 


•  ~TT- 


1  -  A    £Lt 

e^   +  A 


J 


where   A  = 


fllnl 


u1+  flini 


,   u1  =  u(r1),  f11  =  f(r]_,r1) 


These  formulas  are  accurate  if  A   is  snail,  i.e.,  if  the 
magnitude  of  u-   is  greater  than,  or  equal  to,  that  of 


f-j-jn-j,  and  the  accuracy  increases  as 


u 


] 


fllnl 


increases. 


if,  on  the  other  hand,  the  collision  term  is  the 
predominant  one,  we  get  the  following  formulas; 


(8)   n-j  ( t )  =  n,  ee 


-(1  +  v«)   -  u,t 
e 


- ( 1  +  w )   -u2 1 

19)   nr  ( t)  =  n,  ewe         e 


where  u?   =  u(r?),  and  whore  w   is  related  to   t  by  the 
equation 


W 


u 


-ul 


(10) 


1 


nlelll    J 

o 


C      (1  +  w)2,    .    T 

\     q  '    dw  =  1  -  e 


The    integral   ia  tabulated  in   Jahnckc-Emde      (1943     Ed.)    p. 32, 
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v;+l  g2 
since  it  can  be  written  as  /     eS  ds.  (The  use  of  the 


I 


tables   accounts   for   the   seemingly  awkward   form  of   equations 
(8),    (9),    (10)    involving  e.)      These   formulas   are    quite 
accurate    if  u-,    is    of   a    smaller    order    of  magnitude    than 
n-jf-,-,,    and    the    accuracy    increases    as    u-,/f-,-,n-,    decreases. 

As   functions    of  time,    the   foregoing   formulas    are 
especially  accurate   when   t    is    small.      In    section   3,   a    com- 
parison   of    the   formulas   for   a    special   case   with  a   numer- 
ical  solution    illustrates   these    points. 

3 •      Solution  for  a    Special   Case. 

A   solution  has   been   carried   out   in    one    test    case. 
This    special  problem  was    selected    as    a    case    on   the  border- 
line  between   the    two    possibilities   discussed    on    p.    4.      Thus, 
in  any  other  specific   instance,    the  results   will   be    at    least 
as   accurate   as   those    under    consideration  here,    since    the 
accuracy  is    increased   in  either   the   first   or    the    second 
case    according   as    the  ratio    of    the    settling    term    to  the 
collision    term   increases   or   decreases. 

The   example    selected    is        r-,    =  10, n,    =  100. 
Here    u-^   =   110,    f-,,11-,    =  440   so  that    the    collision    term   is 
slightly  predominant   over   the    settling    term. 

This    case   would    correspond    to   the  following   physical 
situation.      Suppose    the  density   P  of    the    liquid    is 
1  gm/cm3   and   the    height   h  is    100   cm.      'The   scaling  factors, 
using    the    units    on   p.    3,    are    then 
<X  =   .083  microns 
1/oc2^     =   1.2  X   106    sec. 


2<*./3/KQ    =  7X108   particles/cm3 
2  oc2^  /K      =  5.8X103   particles/ 
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The    initial   conditions   r-,    =   10,    n-,    =100  mean   an   initial   con- 
centration   of    5.8  X    105  particles/cm3,    all  of   radius 
.85  micron.      The    time  scale    is    such  that    an   interval  of    10" 
in      t      corresponds    to   20  minutes. 

The   equations    (6)    through   (10)   become 

(6)1        nx(t)   =  20/<e110t    -    .8), 


(7)'        n2(t)   -n^O^log 


(10) 


1.394  e^1    -    .394 

•106    '      eDot    +   .894 


(8)'        nx(t)   =  271.8  e 


-(1  +  vO"'   -not 


(9)'        nQ(t)    =  271.8  v;e"(1    +  w)  e"171* 


.0919 


/     e^1   +■* 
do 


'        _  ,  -55t 

dv;   =  1   -  e  * 


We    also   obtained   a    solution   for   this    test    case   by- 
numerical  integration  of    the   differential  equations    (12), 
to  be   derived   in   section  4  from   the   integral  equation    (5), 
neglecting   terms   of   subscript    >  4.      Figure    1  expresses    the    nurn* 
erical   solution    in    graphical  form.      Figure   2   gives    a    pictorial 
comparison   of    the   numerical  solution  for   n^  ( t )    v.ith   the 
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solutions   for  n,   in  (6)',  (8)'»   Figure  3  compares 
the  numerical  solution  with  those  obtained  for   n?   in 
(V)1,  (9)'.   Again  we  repeat,  the  accuracy  is  very  good 
when   t   is  small. 

4«  Mathematical  Treatment* 

In  solving  the  integral  equation  (3)  we  assume,  as 
was  stated  in  section  2,  that  the  initial  distribution  is 
homogeneous.   It  will  be  shown  that  (3)  reduces  to  an  infinite 
set  of  differential  equations,,  of  which  only  the  first  two 
need  be  retained,  the  remaining  ones  being  negligible.   (See 
section  2. ) 

Assume  a  homogeneous  initial  distribution  n(r,0),i.e», 
all  the  particles  have  the  same  radius   r,.   This  can  be 
indicated  mathematically  by  setting 

n(r,0)  =  n1  (j  (r-r^) 

vlicre  6  (r)  is  the   6  -function  and  n,   is  the  initial  concen- 
tration* 

At  any  time   t   the  distribution  will  no  longer  bo 
homogeneous,  but  will  consist  of  particles  of  radius 
r,  ,r?,. . . .  ,r,_,  •  . .  •   Denoting  the  concentration  of  particles 
of  radius   r,   by  n. (t),  the  distribution  n(r,t)   can  then 
be  written  as 

(11)  n(r,t)  =  n^tJo  (r-r1)  +  ng(t)i  (r-r2)  +..+nk(t)  6  (r-rfe)+n  . 

By  substituting  this  In  the  integral  equation  the  following 
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system  of  differential  equations    is   obtained; 


r 


dn1(t) 


oo 


dt 


-n1(t)    [U]_    +    £    fx   n    (t)] 


0=1 


dn2(t) 

dt 


■n0(t)    [u2   +    X  f21n   (t)3    +  |f1]n1(t)n1(t) 


<12)\    dn,(t) 


oo 


ar 

* 

dnk(t) 
— £TE — 


■n, 


-n, 


k 


\ 


t)    [u3   +    Z^n   (t)]    +|f12n1(t)n2(t) 
J-1 


oo 


1=1     _1 


)    [uk+    i:fk.n.(t)]    +  £     g     t^^t^ 


,(t) 


where 


u,     =  u(r7  ) 


and 


fu  =f<ri'rj) 


and   the    initial   conditions    are 


n,(0)    =  n,  ,no(0)   =  n„(0)    =   ...    =  nv(0)    =  0. 


A   closer   analysis   of    the    physical  mechanism  reveals    that 
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particles   of   radius      rn     will   be    immediately  and   rapidly 
created   by   the    coagulation   of    the   original   particles,    but    that 
particles    of  radius      r„      and    larger  are    created    very  slowly. 
As    a   first  approximation,    the   number   of   these   particles    can 
be    neglected;    that    is,    assume 

(13)  n(r,t)   =  n1(t)     I  (r-r-j^  +  n2(  t )     S(r-rvJ. 

Then    (12)   reduces    to 


(14) 


dn,(t) 

— ^  =  -n^t)    [ux   +  f^tt)    +  f12n2(t)    ] 


dn0(t) 


V 


~3T 


=   -n2(t)    [u2   +  f2lnl(t)    +  f22n2(t)]    +  l%Cnl{t^ 


with  initial   conditions   n-,(0)   =  n,,n^,(0)   =  0. 

Before    we    solve    this    system,    it    is   well    to   note    the    fol- 
lowing  relations    between    the    coefficients    u-,  ,u^,f  ^ ,  ,f ,  ?,f  ?p' 
Remembering   that   r2 
can  be    prepared: 


2   

sj   2r-|    =  1.26r-,,    the  following  table 


rl 

X22 

^N 

f12 
fll 

1 

.90 

.97 

1. 

4-2 

2 

.95 

.98 

1. 

48 

4 

.96 

1.00 

1. 

^9 

6 

.97 

1.00 

1. 

54 

8 

.93 

1.00 

1. 

55 

10 

.98 

1.01 

1. 

56 

15 

.99 

1.01 

1. 

57 

oo 

1.00 

1.01 

1. 

59 
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Thus,  for  all  but  the  extremely  minute  particles,  the  coef- 
ficients f,0,f22  c&n  be  safely  replaced  by  f-^,  while  Up  is 
approximately   1.5u-,» 

Case    1,      The    Settling    Term   Predominates. 

In  the   first  equation    (14)   we    can   neglect   f^g^t) 
because   n2(t)    is    small   compared   wiih    n,(t)    in   the   beginning 
and    as      t      increases    it   is    still   small   compared   with  u-,.      This 
leads    to    a   simple   differential  equation   for  n-,(t)   alone,    the 
solution    of   which   is   equation    (6).       (See    p.    5) 

Similarly,    in   the    second   of  equations    (14)    the    term 
fppno(t)    can  be    neglected,    leading    to    the    differential  equation 

nc->     drip(t)  ,  p 

dt        +  n2(t)    [u2    +  f21nl(t)    ]    =  2fllCnl(t,] 

where   formula    (S)    is    to  be    substituted   for    n,(t).      This    is    a 
linear   differential  equation   for  nn(t),    the    solution   of   which 
is 

-/[u9+f91nn(t)]dt     P,  r,     /[uP+fp1n,  (t)]dt 

(16)  n2(t)  =  o  *  .  /^f11[n1(t)]"   e        d      "1J" 

The    quantity  /[uc    +  fQ,n,(t)]dt    can  be    evaluated    immediately 
with    the    result 

/[u2  +f  ni(t)idt     ogt        -»iy2iAn 

e  =  e         L 1   -  Ae  J 

Making   use   of    the    very   close    approximation   fglAll   =  *"» 
formula    (16)  becomes 

-u2t  2   A4         ^   -2U]Lt      u2t 

(17)  n„(t)   = - — -   .  if-..    .   -V-  - ^— - rdt. 


at 


1   -  A2e 


=ujr  •  ^n  ■  -72—     1   - — -p-^r 

11  J  0    1    -   A   e 
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The  integral  can  be  evaluated  explicitly  if  the  further  approx- 
imation u9  =  1.5u-|  is  made.  The  substitution  Ae  ■*■  /  =  x  then 
eivos 


Iv 


iA 


dt 


dx 


1   -  A  e 


9    "ul  L 


UjA 


1   -  x 


u-,A 


log 


1+A 

1    -   Ae-V/* 

1-A    ' 

-u-,t/2 

1    +  Ae      x  '  . 

Inserting   this    in    ( 17 ) ,    one   finally   obtains   equation    (7),    (see 
page    5)   where    the    approximation  u0   =   1.5U-,    has    been    used   for 

-Ur,  t  x 

the    term     e  outside    the    integration   sign. 

Case    2.      The    Collision    Term   Predominates. 

In    this   case    the    preceding  equations  might    hold   at    the 
beginning ,    but  \;ould  be    considerably   in   error   for  large      t. 
It   is    nov;   not    justifiable    to    neglect   f0-,n^(t)    compared   with 

ul  +  fnni{t^# 

Prom   the    table  given    previously,    it    is    seen    that    f .  . 
changes    slightly  when   I   and    j   range   from   1   to   2.      Assume    that 
f..    is    constant   and   represent   it    by  f,« 

The   equations    (14)   are    now 

fdn-,(t) 

— ^f~  =  -n1(t)[u1   +  f1n1(t)    +  f^gtt)] 


(18)< 


dnQ(t) 


Let 


V 


TfF 


-n2(t)[u2   +  f1n1(t)   +  f1n2(t)]    +  ^[n^t )] 


(u-,    -  uQ)t 


n1(t). 


(19)        nn(t)   =  we 

After  some  simplifications,  the  equations  (18)  reduce  to 


/dnn  (t) 


(20)< 


dt 

dw 

dT 


o  (U-,   -  UQ)t 

-u^Ct)  -  f1[n1(t)]^[l  +  v;e  L  d      ] 


1 


f1n1(t)  e 


"(ul  "  u2}t 


V 
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If   tho  first  equation   is   divided   by   the    second,    the 


result   is 


(21)      C^l(t) 

dw 


2u1      (u-.-   u2)t  (u-, 

j —  e  +  2n1(t)[e 


Up)t  2(Un-U0)t 

^      +  we  ]. 


Introduce   the   new    variable 


z   =  e 


(u-l    -   Uo)t 


Then  equation    (21)    and   the    second    of   equations    (20)   become 


(22)  < 


dnx(t) 
cTvv 

dw 
dz 


;u- 


Y~^z    +  2n1(t)    [z   +  viz    ] 


f1n1(t) 


2(u-,    -  u0)z 


2"   • 


If  tho  second  equation  of  (22)  is  solved  for  n-,(t), 
the  result  differentiated  with  respect  to  w  and  then  substi- 
tuted in  the  first  equation,  the  following  equation  for  2 
a   function  of     w      is   obtained: 


as 


(23) 


2    " 

2z    -   Z    Z 


(*') 


M2  u,     -    u 


z   -  £5-    (1   +  wz). 


t 


Primes   denote   differentiation   with  respect    to   w.      The    initial 

conditions   are    z(0)   =  l,z'(0)   =  -  _± __ ,(=   -  k). 

f-1n1iO; 

The    approximation  u9   =  1.5u-,    simplifies    (23)    to 

(24)  z"    =  2z'    z(l   +  wz)    . 

Sot    z    =   1    -   ley  where    y(0)    =  0   and   y'(0)    =   1. 
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Equation  (24)  in  terms  of  y  becomes 

ky"  =  2ky'  (1  -  ky)  [1  +  w(l  -  ky)]. 

If   k  is  small  so  that  it  may  be  neglected  compared  to   1, 
this  equation  reduces  to 

y!1  =  2yf  (1  +  vj)  . 

Its  solution  is 


i  r  (i  +  * ) 

iio 


so  that 


(25)  z  =  1 /  o       '   d\v  . 

°  "° 

Prom  the  second  equation  of  (22),  n,(t)  can  be  found 
to  be 

(26)  nx(t)   =  niea"(1  +  w)  e"ul*  . 

Then  from  equation   (19) 


2    -uot 


c         for     z: 


(27)  n2(t)      =  njowe"*1   +  w )   o 

Finally,    the    relation   between      v»      and      t      is   found   from   (25) 
by  substituting  -, 

loa)    .  "  2*1*    _      ul  r\    (1    +  v;)2     . 

(28)  1   -  e  =   .  a     ■         J     e  avi. 

dilnl       Jo 

In    case    an  expression  for  the    number   n,(t)   of    part- 
icles  of  radius      r„      is    needed,    it   can   be    found   by  solving   this 
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difforential  equation: 


dn3(t) 

dl 


=  -  n3(t)    [u3  +  f31n1(t)    +  f32n2(t)    3    +  f12n1(t)n2(t), 


where  n-,(t),  n2(t)  are  obtained  from  either  (6),  (7)  or  (8), 
(9),  (10).   This  is  a  linear  difforential  equation  which  can 
he  solved  by  standard  methods.   To  evaluate  the  resulting 
integrals,  it  may  be  necessary  to  resort  to  a  graphical  or 
numerical  procedure. 
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